其他
每月一生信流程栏目灵感来自于《铁汉1991》博客的《每日一生信》,他那个时候介绍的主要是生信基础知识,包括数据结构,数据格式,数据库资源,计算机基础等等,所以每天都可以进步,每天都有成果。这些基础知识已经被分享的七七八八了,所以我这里推陈出新,来一个每月一生信流程,陪生信技能树的粉丝们一起进步!
今天学习rnaseqGene
http://www.bioconductor.org/packages/release/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html
在R里面安装这个bioconductor流程
install.packages("BiocManager")
BiocManager::install("rnaseqGene")
1 Introduction 1.1 Experimental data 2 Preparing quantification input to DESeq2 2.1 Transcript quantification and tximport / tximeta 2.2 Quantifying with Salmon 2.3 Reading in data with tximeta 2.4 DESeq2 import functions 2.5 SummarizedExperiment 2.6 Branching point 3 The DESeqDataSet object, sample information and the design formula 3.1 Starting from SummarizedExperiment 3.2 Starting from count matrices 4 Exploratory analysis and visualization 4.1 Pre-filtering the dataset 4.2 The variance stabilizing transformation and the rlog 4.3 Sample distances 4.4 PCA plot 4.5 PCA plot using Generalized PCA 4.6 MDS plot 5 Differential expression analysis 5.1 Running the differential expression pipeline 5.2 Building the results table 5.3 Other comparisons 5.4 Multiple testing 6 Plotting results 6.1 Counts plot 6.2 MA-plot 6.3 Gene clustering 6.4 Independent filtering 6.5 Independent Hypothesis Weighting 7 Annotating and exporting results 7.1 Exporting results 7.2 Plotting fold changes in genomic space 8 Removing hidden batch effects(去除批次效应哦!) 8.1 Using SVA with DESeq2 8.2 Using RUV with DESeq2 9 Time course experiments (时间序列转录组数据)
流程代码
## 表达矩阵来自于R包: airway
# 如果当前工作目录不存在文件:'airway_exprSet.Rdata'
# 就运行下面 if 代码块内容,载入R包airway及其数据集airway后拿到表达矩阵和表型信息
if(!file.exists('airway_exprSet.Rdata')){
library(airway)
data(airway)
exprSet=assay(airway)
group_list=colData(airway)[,3]
save(exprSet,group_list,file = 'airway_exprSet.Rdata')
}
# 大家务必注意这样的代码技巧,多次存储Rdata文件。
# 存储后的Rdata文件,很容易就加载进入R语言工作环境。
load(file = 'airway_exprSet.Rdata')
table(group_list)
# 下面代码是为了检查
if(T){
colnames(exprSet)
pheatmap::pheatmap(cor(exprSet))
group_list
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(exprSet)
# 组内的样本的相似性理论上应该是要高于组间的
# 但是如果使用全部的基因的表达矩阵来计算样本之间的相关性
# 是不能看到组内样本很好的聚集在一起。
pheatmap::pheatmap(cor(exprSet),annotation_col = tmp)
dim(exprSet)
# 所以我这里初步过滤低表达量基因。
exprSet=exprSet[apply(exprSet,1, function(x) sum(x>1) > 5),]
dim(exprSet)
exprSet=log(edgeR::cpm(exprSet)+1)
dim(exprSet)
# 再挑选top500的MAD值基因
exprSet=exprSet[names(sort(apply(exprSet, 1,mad),decreasing = T)[1:500]),]
dim(exprSet)
M=cor(log2(exprSet+1))
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)
# 现在就可以看到,组内样本很好的聚集在一起
# 组内的样本的相似性是要高于组间
pheatmap::pheatmap(M,annotation_col = tmp,filename = 'cor.png')
library(pheatmap)
pheatmap(scale(cor(log2(exprSet+1))))
}
# 以上代码,就是为了检查R包airway及其数据集airway里面的表达矩阵和表型信息
学习这样的流程是需要一定背景知识的
首先是LINUX学习
第1阶段:把linux系统玩得跟Windows或者MacOS那样的桌面操作系统一样顺畅,主要目的就是去可视化,熟悉黑白命令行界面,可以仅仅以键盘交互模式完成常规文件夹及文件管理工作。 第2阶段:做到文本文件的表格化处理,类似于以键盘交互模式完成Excel表格的排序、计数、筛选、去冗余,查找,切割,替换,合并,补齐,熟练掌握awk,sed,grep这文本处理的三驾马车。 第3阶段:元字符,通配符及shell中的各种扩展,从此linux操作不在神秘! 第4阶段:高级目录管理:软硬链接,绝对路径和相对路径,环境变量 第5阶段:任务提交及批处理,脚本编写解放你的双手 第6阶段:软件安装及conda管理,让linux系统实用性放飞自我
然后是R学习
了解常量和变量概念 加减乘除等运算(计算器) 多种数据类型(数值,字符,逻辑,因子) 多种数据结构(向量,矩阵,数组,数据框,列表) 文件读取和写出 简单统计可视化 无限量函数学习
必备书籍及视频
生信技能树关于RNA-seq上下游数据分析的教程的确不少了
转录组经典表达量矩阵下游分析大全 可变剪切,差异转录本,或者差异外显子流程大全 哪怕是到了2018年,RNA-seq仍然可以不做重复 原创10000+生信教程大神给你的RNA实战视频演练 转录组经典表达量矩阵下游分析大全 可变剪切,差异转录本,或者差异外显子流程大全 RNA测序究竟有多可靠呢 TCGA数据库里面的乳腺癌样本RNA-seq数据是配对的有哪些? 生信小白的RNA-seq实战历程 RNA-seq数据分析指南